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Abstract 

We present a systematic computational approach to the study of 
self-similar dynamics. Through the use of what can be thought of as 
a "dynamic pinning condition" self-similarity is factored out, and a 
transformed, non-local evolution equation is obtained. The approach, 
which is capable of treating both first and second kind self-similar 
solutions, yields as a byproduct the self-similarity exponents of the 
original dynamics. We illustrate the procedure through the porous 
medium equation, showing how both the Barenblatt (first kind) and 
the Graveleau (second kind) self-similar solutions naturally arise in 
this framework. We also discuss certain implications of the dynamics 
of the transformed equation; in particular we discuss the discrete-time 
implementation of the approach, and connections with time-stepper 
based methods for the "coarse" integration/bifurcation analysis of mi- 
croscopic simulators. 

1 Introduction 



When one tries to computationally locate constant shape traveling wave 
solutions in an equation that supports them, the standard procedure is to 
go to a traveling frame, so as to change the time-dependent problem into 
a steady state one. Since, however, the correct speed is not, in general, 
known a priori, two possibilities exist: either we will see no steady states 
(having guessed the wrong speed), or we will see infinitely many (if we were 
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fortunate in guessing the correct speed). The simple solution to writing 
fixed point algorithms for this problem is to realize that one is allowed to 
pick one out of the the one-parameter infinity of steady state solutions of 
the PDE in the traveling frame. One-parameter here corresponds to travel- 
ing (translational invariance, a continuous symmetry group) in one spatial 
dimension; the concept generalizes directly to traveling in more than one 
spatial dimensions. Computationally one augments the system by writing 
an additional (possibly nonlocal) condition, called a "pinning condition" [^. 
We then have enough equations to solve for the one solution and the correct 
traveling speed simultaneously. 

Rowley and Marsden [10| introduced recently a template-based technique 
for "factoring out" translational invariance; their work was motivated by the 
so-called Karhunen-Loeve expansion in the case of PDEs with traveling wave 
solutions. This was done by choosing a (more or less arbitrary) function, 
the template, and using it to systematically pick out one out of the one- 
parameter infinity of elements of the group orbit in a manner that will be 
explained below. We have recently implemented a discrete time version of 
this template-based approach to perform what we call "effective" bifurcation 
analysis of evolution equations in complicated media [11 1. 

An important contribution of the work of Rowley and Marsden is that 
through the template method they allow the pinning to be done in a dy- 
namic way. This means that we do not just implement some computational 
contraction mapping, and then converge on its fixed point, which is the 
solution traveling with constant speed and constant shape. Through the 
template we actually have a way to "factor out" the translational invariance 
in continuous time, and thus get an equation in an appropriate frame at all 
times. This is quite interesting since, often, what we want to find is not just 
the final solution, but we also want a view -"unobstructed" by the symme- 
try group- of the physical dynamics of the approach to the traveling wave. 
It is of course, in principle, possible to first find the traveling solution by 
some independent means (e.g. analytically) and then factor it out from the 
dynamics. The template approach allows us to do the "factoring out" of the 
group orbit naturally and continuously in time, without a priori knowledge 
of either the traveling solution or its speed. Moreover, as a simple byproduct 
it yields at the end the solution and the correct speed. 

In this paper we extend the template approach to another continuous 
symmetry group: self-similarity. We take advantage of the fact that using a 
template to factor out self-similarity, we actually obtain a dynamic evolution 
equation. This equation not only gives the self-similar solutions, but it 
actually gives a view of the dynamics of approach to the self-similar solution 
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with the self-similarity conveniently factored out. 

The remainder of the paper is organized as follows: in the following 
section we show how to derive the transformed evolution equation for the 
scaled variable. We then illustrate the procedure through the computation 
of two well-known solutions (the Barenblatt solution, a self-similar solu- 
tion of the first kind, and the Graveleau solution, a self-similar solution of 
the second kind) through the dynamics of the transformed equation. We 
conclude by a discussion of what we perceive as the implications of the ap- 
proach, including alternative (to integration) computational approaches to 
the location of self-similar solutions. We also discuss briefly the discrete- 
time implementation of the approach on the original (as opposed to the 
transformed) equation. This has possibly important implications because of 
the recently proposed "coarse" integration/bifurcation analysis techniques 
for analyzing microscopic (e.g. Molecular Dynamics (MD) or Monte Carlo 
(MC)) timesteppers. It is conceivable that the computer-assisted analysis of 
self-similarity proposed here can be carried out even in cases where coarse, 
macroscopic equations describing the evolution of moments of molecular 
distributions conceptually exist but are not available in closed form. 

We also note that what we propose here shares many elements with the 
approaches described, for example, in page 326 of the book of Goldenfeld 
, and in Chapter 6 of , including references to the work of the authors 
of the book, as well as Papanicolaou, Zakharov and their groups. We are 
currently exploring the similarities and differences of our work with that 
literature (see also Q). 

Because of the long ties, over the years, of two of the authors with the 
University of Minnesota, and because this work came to fruition in an office 
at U of M., we will call the transformed dynamics developed here the "MN- 
dynamics" . 

2 Self-similar solutions 

Self-similar solutions play an important role in the development of the the- 
ory of nonlinear evolutions equations. In addition to providing exact and 
sometimes even explicit solutions they often describe the asymptotic behav- 
ior at large times, or the asymptotic form of solutions in the neighborhood 
of some important change in behavior such as focusing or blow-up at some 
finite critical time. Often a self-similar solution to an autonomous evolution 
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equation in the variables {xi,X2, ...,Xn,t) is a function of the form 

It - t*r' ' It - 1*\ 



\t-t*\PF — — (1) 



where the similarity exponents ai, and the function F must be de- 

termined from the equation together with appropriate boundary and initial 
conditions. Here t* is either the known initial time in which case we con- 
sider t > t* , or the unknown critical time so that we consider t < t* . In 
view of the time translation invariance we can always take t* = 0. In some 
problems the exponents can be obtained a priori from scaling arguments 
and conservation laws. This is referred to as self-similarity of the first kind. 
However, it is often the case that the exponents cannot be gotten a priori 
and are usually obtained by solving what amounts to a nonlinear eigenvalue 
problem for the function F. This is the case of self-similarity of the second 
kind. A cogent account of the theory of self-similar solutions with many 
illuminating examples can be found in Barenblatt's book Q. 

In the next section we will derive the MN-dynamics equations for the 
construction of self-similar solutions. In section 4 we illustrate the method 
by applying it to the computation of the Barenblatt solution (first kind self- 
similarity) and the Graveleau solution (second kind self-similarity) of the 
porous medium equation (pme) 

dt ^ ' 

where m > 1 is a constant. For properties of the porous medium equation 
see in. 



3 MN-Dynamics 

Consider the partial differential equation 

ut = iu) , (2) 

where the generally nonlinear operator D acts on the variable x. We restrict 
ourselves to operators for which there exist two constants a and h such that 
for all positive A and B 

(Bf (^^)) = A-B'Dy if (y)) , where y = j. (3) 
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For these operators the PDE may have a one-parameter family of self-similar 
solutions 

uix,t) = s^U i^-^y (4) 
where C > is the parameter, s = \t — t*\, and U satisfies the ODE 

ctC~" {pU - ayUy) = Dy {U) in the variable y = (5) 

with a = sgn(t — t*). The exponents a and f3 satisfy the scaling condition 

p-l = aa + bp. (6) 

We now present an alternative method of computing similarity solutions, 
which factors out the similarity exponent a. We start with the general 
scaling 

u{x,t) = Bis)w(^-^yT{s)y (7) 
where A,B, and r are unknown functions. Using Eq. (H) we obtain the 



PDE for w 



^ ( - ^-jywy + TsWr ] = A^B^-^Dyiw) (8) 



Given an appropriate template function T{y) we want to adjust A{s) so that 
w satisfies the orthogonality condition 

wiy,T)T{y)dy = 0. (9) 

-oo 

Multiplying Eq. Q by T{y) and integrating we obtain 

^ = -aA'^B'-'lir), (10) 

where 

j^^^ ^ I^^Dyiw{y,r))T{y)dy ^ ^^^^ 
J^^yWy{y,T)T{y)dy 

To determine B we need an additional pinning condition. For numerical 
convenience we impose the local condition w{p,t) = 1 at some fixed point 
X = p. For example, p may be a symmetry point where Wy{p, r) = for all 
r. Evaluating Eq. (^) at y = p and substituting from Eq. ( [ToD we obtain 

^ = aA^B'-' {Dy{w{p, r)) - pwy{p, r)/(r)) . (12) 
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Substituting Eqs. (|ToD and in Eq. (|) we obtain 

Wr + /(r) {ywy - wpwy{p, r)) + wDy{w{p, r)) = Dy{w) (13) 
provided that 

Ts{s) = aA''{s)B'-\s). (14) 

In view of Eq. ([l^ we can rewrite A and i? as functions of r instead of 
s. Let A{t) = A{s{t)) and B{t) = B{s{t)). Then it follows from Eqs. 
(ID, and dH) that 

= -H^) and ^ = Dy{w{p,T)) - pWy{p,T)I{T). (15) 

Moreover, the physical time s = |t — t*| is given as a function of r by 

^ = aA-''{T)B^-\T). (16) 
ar 



Starting from a suitable initial condition w{y,0), we integrate Eq. ( [T^ ) for 
T ^ oo. If w{y,T) converges to a steady state Vl^(y), there is a similarity 
solution to Eq. (^). In particular, W satisfies 

k (yWy - pWyip)W) + WDy(Wip)) = Dy(W), (17) 

where 

k = lim I(t) (18) 

T — ^OO 

as well as Eq. (|5|). Comparing coefficients in Eqs. (^) and ( [T7| ) we find that 

a — 



/3 I)j;I^(p) - kpWyip) 



(19) 



Thus the similarity exponents a and /3 are determined by the linear system 
consisting of Eqs. @ and (|l9|). 

In problems where there is a finite critical time t* we are interested in 
t < t* so that a = —1 and s = t* — t. We integrate Eqs. (p^) and (|^) to 
obtain 

t* = t(T = 0) + lim r A-''(rj)B^-\7])dri. (20) 

A remarkable property of this method is that is able to capture both 
self-similar solutions of the first and second kind, as we show in the next 
section. 

We note that the new PDE Eq. ( [l3|) can also be obtained by repeatedly 
evolving Eq. (^) an infinitesimal time interval dt, and then rescaling the x 
and u variables such that they satisfy u{L,t) = 1 and Eq. (^). 
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4 Examples 



Here we give two examples of numerical integration of Eq. (13) to obtain 
self-similar solutions to the porous medium equation with m = 2. 

The ID Barenblatt solution is the self-similar solution of the first kind 
to the porous medium equation whose initial datum is a Dirac mass at the 
origin. This solution is known explicitly and 

a = -(3= 1/3. 

The operator in this case is Dx{u) = {v?')xx and the scaling condition Eq. 
(^) reduces to /? = 2a — 1. We discretized the equation with centered 
finite differences. The time integration is carried out with an explicit Euler 
scheme on a domain [0, 3] discretized with 600 points. From the linear 
stability analysis of the numerical scheme, the variable time step must satisfy 
At < 0.5/|w(y, t)| on the domain of integration. To compute the Barenblatt 
solution we choose the template function 



ny) 



1 for \y\ < 1 
— 1 for \y\ > 1 



w{y,0) 



The orthogonality condition Eq. (^) means that the dimensionless mass 
within \y\ < 1 equals the mass outside this interval. For the initial datum 
we take 

1 for \y\ < 2 
for \y\ > 2 

which satisfies the orthogonality condition. Fig. 1 shows the convergence 
of the numerical solution to the steady state. The numerical value of the 
similarity exponent given by Eq. ( [T^ ) is a = 0.333371 for r = 10. 

The Graveleau solutions to the porous medium equation form a 1-parameter 
family of axi-symmetric focusing solutions which are self-similar of the sec- 
ond kind Q. To construct them by integrating the MN-dynamic equation 
we write the equation in (r, t) coordinates. For axi-symmetric diffusion the 
operator is Dr{u) = ( 

vP^rr ~l~ (^^)r/^ and the scaling condition is again 
P = 2a — 1. We approximate the solution that has an infinite support with 
a numerical solution in a finite interval [0,10], discretized with 2000 grid- 
points. We use the boundary conditions 'w{0,t) = 0,Wy{L,T) = 0, and in 
Eq. (|1|) we set L = 10. 

The template function is chosen to be 



\ 1 for 7 < \y\ < 10 
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and the initial datum is 

, f for < lyl < 4 

^^^'°) = \ 1 for4<|y|<10 

which is orthogonal to T{y) on the integration domain. We show the evolu- 
tion of w in Fig. 2. As r ^ oo the solution tends to a steady state. In this 
case the corresponding physical time tends to a finite value t* which is the 
time at which the solution to the original initial value problem first becomes 
positive in the neighborhood of the origin, t* is called the focusing time. 
To determine the similarity exponent we use Eqs. (^) and (19). For 2 space 
dimensions the correct value of a = 0.856333... and the MN-dynamics gives 
the approximate a = 0.85695. The accuracy would be improved by using 
more gridpoints and enlarging the domain. 



5 Discussion 



Locating the self-similar solutions and their exponents by integration of the 
MN-equation is, of course useful. It is, however, only one of the ways that we 
can now find these solutions. Now that the self-similar solutions have become 
steady states of the MN-equation, we can bring all the tools of nonlinear 
dynamics to bear on their study. In particular, we can find self-similar solu- 
tions not only through the fact that they may be attracting (by integration), 
but also through contraction mappings like Newton's method (this would 
be interesting in physical problems like the Landau-Gizburg equation fisll ). 
Indeed, we can do bifurcation / continuation / stability calculations of the 
closed non-local PDF and find all types of stable, unstable, bifurcating, pos- 
sibly even more complicated (e.g. periodic in the transformed frame, or even 
more complicated, chaotic in the transformed frame) self-similar solutions. 
They are going to be elements of the global attractor of the MN-equation, 
or, in general, invariant objects for the MN-dynamics. For example, in the 
traveling limit cycle in the transformed equation would correspond 

to a modulated traveling wave — so, possibly, a limit cycle of the trans- 
formed equation would correspond to a "periodically self-similar" solution 
for the self-similar problem. These might arise in Hopf bifurcations of the 
MN-equation, and might be associated with complex exponents; we are not 
sure at the moment what all these additional invariant objects might be, 
but it will be interesting to study them. We are certainly aware of cases 
that should correspond to steady state bifurcations in the MN-dynamics. 

More importantly, however, the Jacobian of the linearization of the MN- 
equation at its steady states or limit cycles, should it prove to be a well 
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defined object, will tell us what the stability in the self-similar frame is, with 
the self-similarity factored out. It will help us understand asymptotic rates 
of approach, bifurcations of new self-similar solutions, possibly find basins 
of attraction of different self-similar solutions as stable manifolds of saddle- 
self-similar solutions. In particular, it is conceivable that in the case of 
Type-II self-similar solutions, continuation and the use of this lenearization 
might facilitate the study of the post-focussing regime. There are many 
interesting singularities whose "other side" we would like to explore. The 
porous medium equation is one of them; possible examples from cosmology 
[13 1 also come to mind. 

We conclude with a few (speculative at the moment) comments. 

First, when a successful transformation turns the problem of finding self- 
similar solutions into finding steady states, limit cycles, chaotic attractors, 
or in general, elements of the global attr actor of the transformed dynamics, 
we can bring all the machinery of dynamical systems / bifurcation theory to 
bear on computing self-similar solutions and their transitions (bifurcations 
in the tranformed frame). In particular, techniques much more powerful 
than simple integration (e.g. contraction mappings like the Newton method, 
continuation, and multiparameter bifurcation theory) can be brought to bear 
on the transformed problem. 

Second, in everything we have done so far, we first transform the equa- 
tion, and then work with the transformed equation. In effect, one has to 
write a new code to analyze the old problem. The original reason we were 
attracted to the template approach to factoring out translational invariance 
was that we were interested in using it in what in general we term "time- 
stepper" based methods for numerical bifurcation theory. It is possible to 
implement the template procedure not only in continuous time, but also 
in discrete time: one runs the "original" dynamics for some time and then 
"pulls back" the result (as opposed to constantly pulling it back as in |Q]). 
We have implemented this discrete time, time-stepper based bifurcation ap- 
proach in [11 1 for the traveling case; in the same spirit, the discrete-time 
version of the approach can be used on the original equation in the self- 
similar case (as opposed to writing down the transformed MN-equation). 

This opens the way for another interesting possibility. Over the last 
few years, we have developed methods for what we call "coarse bifurcation 
analysis" of microscopic time-steppers [12, |^, ^. In these, we create a map 
(using microscopic dynamics, such as Monte Carlo, or Lattice Boltzmann) 
from an initial condition in the space of a few moments of a distribution to 
a final condition (after some time) in the same truncated moments space 
of the distribution. The procedure involves "lifting" the macroscopic initial 
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condition to one or more microscopic ones conditioned on the few governing 
moments, running the microscopic code, and then averaging (restricting) 
back to governing moments space. There is a number of interesting issues 
about how long to run the microscopic simulator, and how to do good vari- 
ance reduction, but we will not discuss this here. 

Based now on the coarse-time stepper implementation, suppose that we 
have microscopic simulations (simulations at one level of description), and 
we suspect that some coarser description (moments) of the microscopic dy- 
namics will have self-similar behavior. Then, using the coarse time-stepper 
and the discrete-time version of the above procedure, it is possible, given a 
set of reasonable assumptions discussed in some detail in [0] to find cer- 
tain elements of the self-similar dynamics (e.g. steady states, invariant sets, 
bifurcations) of what the "coarse-MN-equation" would have been, without 
having to explicitly construct (approximate) this coarse-MN equation. In- 
deed, this approach gives rise to alternative, mathematics-motivated ensem- 
bles for performing the microscopic simulations (with an eye towards the 
analysis of coarse dynamics) |^]. Recent advances in multiscale computa- 
tions (projective and telescopic projective integrators (^,0, as well as what 
we provisionally call "patch dynamics") may assist in doing these calcula- 
tions efficiently, in smaller space and time scales than the full computational 
domain. These developments in multiscale computations are rather general, 
and not designed in particular with the coarse MN-dynamics in mind. 

We believe that the general approach outlined in this paper opens new 
computational possibilities in the study of self-similar problems and their 
dependence on parameters. Problems ranging from interfacial fluid dynam- 
ics to shock waves, and from cell aggregation to cosmology and materials 
science may become more accessible to computer-assisted analysis, whether 
through closed macroscopic equations, or through the "coarse" analysis of 
alternative, more microscopic descriptions. 
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Initial condition 



0.55 




-0.55 



Figure 1: Evolution of w{y,T) for the porous medium equation. The ini- 
tial rectangular profile evolves towards the steady state, keeping the scalar 
product with the template (dashed line) equal to zero. 
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Figure 2: Evolution of w{y, r) for a convergent flow in the 2D porous medium 
equation. The initial rectangular profile evolves towards the steady state, 
which describes the selfsimilar solution before the focusing time. 
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